#include "nr.h"

void calc_NR(int N_NODES,
             int N_ELEMENTS,
             int N_ELECTRODES,
             double HEIGHT_LAYER,
             double **COORDS,
             int **TOPOLOGY,
             double *V_RHO_0,
             int CURR_PATTERN,
             int *V_ELECTRODES,
             double **VOLT_PATTERN_MATRIX,
             double **REGULARIZATIONS_MATRIX,
             int  *RtR_SIMPLIFY,
             double VALOR_CORR,
             int no_terra,
             int N_NODES_LOCAL,
             int *VARIABLE_TABLE,
             double *volt_medido,
             int FLAG_MODELO_ELETRODO_3D,
             int *erro_saida)
{
  double **k_global,***k_locals,**H;
  double *volt_calculado;

  volt_calculado=dvector(1,N_ELECTRODES*N_ELECTRODES);

  switch (N_NODES_LOCAL){
    case 3:
      monta_rigidez_tri(&k_global,
                         &k_locals,
                         N_NODES,
                         N_ELEMENTS,
                         N_ELECTRODES,
                         HEIGHT_LAYER,
                         COORDS,
                         TOPOLOGY,
                         V_RHO_0,
                         N_NODES_LOCAL,
                         VARIABLE_TABLE,
                         erro_saida,
                         FLAG_MODELO_ELETRODO_3D);
      break;
    case 4:
      monta_rigidez_tetra(&k_global,
                           &k_locals,
                           N_NODES,
                           N_ELEMENTS,
                           N_ELECTRODES,
                           COORDS,
                           TOPOLOGY,
                           V_RHO_0,
                           N_NODES_LOCAL,
                           VARIABLE_TABLE,
                           erro_saida,
                           FLAG_MODELO_ELETRODO_3D);
      break;
    case 6:
      monta_rigidez_penta(&k_global,
                           &k_locals,
                           N_NODES,
                           N_ELEMENTS,
                           N_ELECTRODES,
                           COORDS,
                           TOPOLOGY,
                           V_RHO_0,
                           N_NODES_LOCAL,
                           VARIABLE_TABLE,
                           erro_saida,
                           FLAG_MODELO_ELETRODO_3D);
      break;
    default:
      (*erro_saida)=302;
      break;
  }

  if((*erro_saida)!=0){
    printf("\nERRO numero %d - %s.\n",*erro_saida,"erro no monta_rigidez da fase 1");
    free_dmatrix(k_global,1,N_NODES,1,N_NODES);
    free_d3tensor(k_locals,1,N_ELEMENTS,1,N_NODES_LOCAL+1,1,N_NODES_LOCAL);
    free_dvector(volt_calculado,1,N_ELECTRODES*N_ELECTRODES);
    return;
  }

   calcula_jacobiano(N_ELEMENTS,
                    N_ELECTRODES,
                    N_NODES,
                    CURR_PATTERN,
                    VOLT_PATTERN_MATRIX,
                    VALOR_CORR,
                    no_terra,
                    k_global,
                    k_locals,
                    V_ELECTRODES,
                    NULL,
                    N_NODES_LOCAL,
                    V_RHO_0,
                    &H,
                    volt_calculado,
                    VARIABLE_TABLE,
                    erro_saida);

  if((*erro_saida)!=0){
    free_dmatrix(k_global,1,N_NODES,1,N_NODES);
    free_d3tensor(k_locals,1,N_ELEMENTS,1,N_NODES_LOCAL+1,1,N_NODES_LOCAL);
    free_dvector(volt_calculado,1,N_ELECTRODES*N_ELECTRODES);
    free_dmatrix(H,1,N_ELECTRODES*N_ELECTRODES,1,N_ELEMENTS);                VARIABLE_TABLE[20]=0;
    printf("\nERRO numero %d - %s.\n",*erro_saida,"erro no calcula_jacobiano");
    return;
  }

  atualiza_rho (N_ELEMENTS,
                N_ELECTRODES,
                H,
                V_RHO_0,
                REGULARIZATIONS_MATRIX,
                volt_calculado,
                volt_medido,
                RtR_SIMPLIFY,
                erro_saida);

  if((*erro_saida)!=0){
    free_dmatrix(k_global,1,N_NODES,1,N_NODES);
    free_d3tensor(k_locals,1,N_ELEMENTS,1,N_NODES_LOCAL+1,1,N_NODES_LOCAL);
    free_dmatrix(H,1,N_ELECTRODES*N_ELECTRODES,1,N_ELEMENTS);                VARIABLE_TABLE[20]=0;
    free_dvector(volt_calculado,1,N_ELECTRODES*N_ELECTRODES);
    printf("\nERRO numero %d - %s.\n",*erro_saida,"erro no NR_fase1");
    return;
  }

  free_dmatrix(k_global,1,N_NODES,1,N_NODES);
  free_d3tensor(k_locals,1,N_ELEMENTS,1,N_NODES_LOCAL+1,1,N_NODES_LOCAL);
  free_dmatrix(H,1,N_ELECTRODES*N_ELECTRODES,1,N_ELEMENTS);                VARIABLE_TABLE[20]=0;
  free_dvector(volt_calculado,1,N_ELECTRODES*N_ELECTRODES);
}
